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Abstract 

A thermal lattice Boltzmann model is constructed on the basis of the ellipsoidal statistical 
Bhatnagar-Gross-Krook (ES-BGK) collision operator via the Hermite moment representa- 
tion. The resulting lattice ES-BGK model uses a single distribution function and features 
an adjustable Prandtl number. Numerical simulations show that using a moderate discrete 
velocity set, this model can accurately recover steady and transient solutions of the ES-BGK 
equation in the slip-flow and early transition regimes in the small Mach number limit that 
is typical of microscale problems of practical interest. In the transition regime in particular, 
comparisons with numerical solutions of the ES-BGK model, direct Monte Carlo and low- 
variance deviational Monte Carlo simulations show good accuracy for values of the Knudsen 
number up to approximately 0.5. On the other hand, highly non-equilibrium phenomena 
characterized by high Mach numbers, such as viscous heating and force-driven Poiseuille 
flow for large values of the driving force, are more difficult to capture quantitatively in the 
transition regime using discretizations chosen with computational efliciency in mind such 
as the one used here, although improved accuracy is observed as the number of discrete 
velocities is increased. 



1 Introduction 



The lattice Boltzmann (LB) method has been successful in simulating isothermal (athermal) 
flows; however, thermal flows still remain a challenge, despite significant efforts from a number 



of research groups (see 


Shim & Gatignol 2011; Sbragaglia et al., 2009j |Scagliarini et al. 


|2010 


Prasianakis et al. 2009 


Prasianakis & Karlin, 2008j |Watari, ,2009; jConnella et a/.||2007|l 


Zheng 


et al. 


2010 Guo et al. 


20071 Shan & Chen, 2007| Watari & Tsutahara, 2003| Prasianakis & 


Karlin 


2007 Lallemand & Luo 2003 Chen & Doolen 1998). Recent approaches fall into two 


broad categories: the high-order approach (e.g., Sbragaglia et al. 2009| Shan & Chen 


2007 


Watari & Tsutahara 20031 and the double-distribution- function approach (e.g.. He et al. 


1998 



Shan 1997). As a direct extension of the isothermal LB model, the high-order model only 
uses one distribution function while retaining high-order terms in the equilibrium distribution 
function to recover the full Navier-Stokes (NS) equations. Unfortunately, this approach requires 
a richer velocity space and as a result suffers in terms of simplicity and in some cases numerical 
instabilit}|^ By contrast, in the double-distribution-function approach, two distribution functions 
are utilized: one for the velocity field and the other for the temperature field. As a result, high- 
order terms can be avoided and the standard lattice can be used. In addition, this approach is 
numerically more stable. However, to utilize the standard lattice, the energy and momentum 



^Wo also note some efforts on constructing single-distribution-function thermal models using standard lattices 
e.g. a 2 dimensional model with nine discrete velocities (see [Prasianakis fc KarlTn] |2007"| l . 



transport equation need to be decoupled through the Boussinesq approximation (cf. Guo et al. 



2007); in this case the double distribution approach does not recover the full NS equations in a 



strict sense. 

In addition to heat transfer, kinetic effects can make flows even more complex. For example, 
for gaseous flows at the microscale, kinetic effects have to be taken into account as the Knudsen 



number (Kn, the ratio of the mean free path and the characteristic length) becomes finite ( Had- 



jiconstantinou 2006). Due to its kinetic origins, the LB method has been shown to be able to 



describe moderately complex kinetic effects (see, e.g., Zhang et aZ.', '2005 Toschi & Succil 



et al, 2007; Ki m et al\ 2008 a| Yudistiawan et al. 



2005 



Sbragagha k Succi[,200 5, 2006 ; Tang et aL[|20086| [Zhang et al, 2006; Shan aLH2006, Ansumah 



20081 |Guo erarr|2006||Succi||2002| |Kim et al 



20086; |Tang et al 2008a Tian et al 2007). Particularly, as the LB method can be considered 



as an approximation to the Boltzmann-BGK equation, high-order models should, in principle, be 
able to go beyond NS hydrodynamics ( Shan et al 2006 Meng & Zhang 2011a ). Both analytical 



and numerical analysis have shown that using an appropriate discrete velocity set, LB models can 
capture kinetic effects including the Knudsen layer qualitatively and provide reasonably accurate 
results for a range of isothermal problems (Kim et al 2008a Yudistiawan et al 2008); also, 
as expected, more discrete velocities generally give better predictions ( Meng k, Zhang, ,2011 a| . 
When gas-surface interactions are concerned, discrete velocity sets from even-order Hermite poly- 
nomials typically perform better (Meng & Zhang 20116). For thermal rarefied flows, despite 
that the velocity-slip and temperature-jump problems have been investigated using a high-order 
LB model (Watari 2009), signiflcant effort is still required to develop robust LB methods for 



modeling general thermal non-equilibrium flows of engineering interest. 

The Hermite expansion provides a systematic approach to derivation of the LB models for 
representing the NS hydrodynamic systems and beyond (Shan et al 2006). The usage of the 



Hermite expansion was pioneered by Grad ( 1958 1949 ) for approximating solutions to the Boltz- 
mann equation. An important feature of Hermite polynomials is that the expansion coefficients 
correspond directly to moments of the distribution function. In addition, the truncation of 
higher-order terms does not directly alter the lower-order velocity moments of the distribution 
function. In his famous paper, Grad kept the first thirteen Hermite coefficients and obtained 
the well known 13-moment system (Grad 1949). This system contains physics beyond the NS 



equation. However, a significant weakness of Grad's 13-moment system is its inability to describe 
smooth shock structures above a critical Mach number. Also, this system is not easily related 
a priori to the Knudsen number. As a result, significant effort has been devoted towards the 
development of an improved moment system in the last decade. Regularization for Grad's 13 
moment system (i.e., the so-called R13 system) has been introduced on the basis of a Chapman- 
Enskog expansion around a non-equilibrium state (Struchtrup & Torrilhon 2003 Struchtrup 



2005 2004 ) . The related issues associated with the numerical scheme and boundary conditions 
have then been discussed bylTorrilhon & Struchtrupl (|2008p andlGu & Emersonl (|2007p . A series 



of analytical solutions were also obtained and compared to the kinetic solutions (see Taheri et al 



2009 Taheri & Struchtrup 2010). In addition to the R13 system, higher-order moment equation 



models have also been developed ( Gu & Emerson 2009 ) . 

Inspired by Grad's work, the Hermite expansion has been used to derive new LB models 

It was shown that the truncation of the Hermite expansion 



Shan et al 


2006 


Shan & He 


1998 



is equivalent to evaluating the distribution function at the chosen order of the discrete velocities. 
Thus, the resulting LB models can describe gaseous systems at the corresponding level. However, 
contrary to Grad's approach, the governing equation of LB models is presented in a much simpler 
form (e.g., the linear convective term) as the evolution is accomplished at the distribution function 
level rather than the state- variable level. Therefore, the LB model has some advantages, e.g., 
straightforward numerical implementation and code parallelisation. 

In this work, we establish a thermal LB model based on the ellipsoidal statistical Bhatnagar- 
Gross-Krook (ES-BGK) equation following the systematic procedure of high-order Hermite ex- 
pansion. The resulting model features an adjustable Prandtl number in contrast to the BGK 
model whose Prandtl number is fixed at 1. We validate our model using shear-driven (Couette) 
flows, Fourier flows and unsteady boundary heating problems in the small Mach-number limit 
over a range of Knudsen numbers. Using only a moderate set of discrete velocities, we flnd 
very good agreement with benchmark solutions. We also show that, with the same moderate 
discrete velocity set, highly non-equilibrium problems characterized by high Mach numbers, such 
as force-driven Poiseuille flow for large values of the driving force, are captured accurately in the 



2 



slip-flow regime and with good qualitative accuracy in the transition regime. Better accuracy 
requires a higher number of discrete velocities, in addition to higher-order Hermite expansion. 



2 Thermal lattice Boltzmann model 
2.1 The ES-BGK equation 

The difficulty associated with solving the Boltzmann equation is mainly due to the collision 
term. To reduce the complexity, the collision term may be replaced by a simple statistical model 
that preserves the conservation laws of mass, momentum and energy. The most commonly-used 
model is known as the Bhatnagar-Gross-Krook (BGK) collision operator, which is often used 
in constructing lattice Boltzmann models. Despite its success, the BGK model suffers from a 
number of drawbacks: specifically, it yields a fixed Prandtl number (Pr) of unity, while the 
correct value for a monatomic gas is close to 2/3. To address this limitation, Holway' f 1966|) 
proposed the ES-BGK model which replaces the Maxwellian distribution of the standard BGK 
equation with an anisotropic Gaussian distribution. This model can be written as 



dl 
di 



a/ 

'' dxi 



9i 



do, 



1 



(1) 



where / denotes the single-particle distribution function, the phase (particle) velocity, cji the 
body force and f the mean relaxation time which is assumed to be independent of particle 
velocity. The anisotropic Gaussian distribution can be written as 



fES 



1 



det[27rAi 



: exp 



1 



(2) 



where Xij = RTSij + {haij)/ p and R is the gas constant. Macroscopic quantities, such as density 
p, velocity it, shear stress ct, temperature T and heat flux q, can be obtained by integrating the 
distribution function over the velocity space, i.e. 



p 
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pUi 




Ci 




-!> 




ft 




1 

2 i J J 


pDRT 




Ci Ci 



dc. 



(3) 



where Ci = ct — Ui is the particle peculiar velocity and the angle bracket denotes the trace-free 
tensors (see Struchtrup 20056 Appendix A. 2. 2). The pressure p can be related to the density 



and temperature by the equation of state 



p = pRT. 



(4) 



^6^1. The 



As (the inverse of the matrix) must be positive definite, h is restricted to — \ 
exact values of b and f can be determined by the Boltzmann integral with a known inter-molecular 
force law or experimental data. In this work, experimental values will be used. The Navier Stokes 
equations can be derived (Holway 1966) using a Chapman-Enskog expansion. The viscosity 
and thermal conduction coefficients can be written as jl — {pf)/{l — h) and k — pR{D + 2)f/2, 
respectively, where D is the spatial dimension. Therefore, the Prandtl number is given by 
Pr — 1/(1 — b) and can be adjusted via the free parameter b. For thermal flows, viscosity 
depends on temperature, and can be expressed as fi/fio = (T/Tq)^, where w is related to the 
molecular interaction model, varying from 0.5 for hard-sphere molecular interactions to 1 for 
Maxwellian interactions. Therefore, in general, f may be taken to be a function of temperature. 

The ES-BGK model has received renewed attention, in part because the associated H theorem 
was recently proven (Andries et al. 2000). Several studies have shown that it can provide 



reasonable accuracy for modeling transport in simple rarefied flows (Andries et al. 



2002 



Han 



eTaTl |2007| |Graur fc Polikarpov| [20091 |Han etlU^Un] [Mieussens fc Struchtrupj |2004[ IGalUs 



& Torczynski 2011 Garzo &: Andres Santos 1995). In the following, we will develop a thermal 



LB model based on the ES-BGK equation using the Hermite moment representation (Shan 
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et al. 2006); we will refer to the resulting LB method as the lattice ES-BGK model. In order to 



separate the numerical error arising from the LB method from the modeling error associated with 
approximating the hard-sphere operator with the ES-BGK model, we will also present numerical 
solutions of the Boltzmann equation with the ES-BGK collision operator obtained using the 



finite difference/discrete velocity method presented in Aoki et al. (2002 1. Although this method 
is not practical for problems of engineering interest due to the high dimensionality associated 
with the very fine discretization of the velocity space, it can be used to obtain accurate solutions 
of the ES-BGK equation in the one-dimensional benchmark problems considered here. In order 
to differentiate from our lattice-based solutions, we will refer to the ones obtained by the finite 
difference/discrete velocity method as numerical. 

2.2 Derivation of the lattice ES-BGK equation 



Shan et al. ( 2006 1 have shown that a hierarchy of different order LB models can be systematically 



derived from the BGK equation via the Hermite expansion, which can be regarded as approxi- 
mations to the BGK equation. Gas flows can then be described at the kinetic level with various 
discrete velocity sets. This procedure can also be used to derive the lattice ES-BGK model. For 
convenience, we use the following non-dimensional system. 



L 



,t 



VRT oi 
L 



Lgi 



1 



L Po Po Po Po PovRTo 

The temperature and state equations become 

DT =- 
P . 

P = pT. 

The non-dimensional relaxation time can be written explicitly as 

poy/RTo p Kn p 



T 
To' 



Po 



fCiCidc, 



PqL Prp Pr p ' 



where the Knudsen number is defined by 



Kn 



PoVRTo 
PoL 



To solve Eq.([T]), the distribution function is first projected onto a functional space spanned 
by the Hermite basis. If a Gauss-Hermite quadrature of a degree ^ 2A^ is chosen, then the 



first A'' velocity moments of the distribution function are retained (see Shan et al. 2006 P.420). 
Therefore, the distribution function is essentially approximated as 



^ 1 

fix, c, t) « Lo{c) J2 ^a^"^ i^, t) ■■ X^"Hc), 

n=0 

and the coefficients a'"^ are calculated from 

where x^"^ is the n-th order Hermite polynomial and w(c) is the weight function. 
The key step is to expand the anisotropic Gaussian distribution as 

71=0 



(5) 



(6) 



(7) 
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where 



41 = / fESX^-^dc^Y^ ^/|'^^"(c„), 
a=l ^^^°') 



(8) 



and Wa and Cq, a = 1, • • • , d, are the weights and abscissae of a Gauss- Hcrmitc quadrature. For 
the anisotropic distribution, the first a few coefficients are given by 



a 



(3) 

ES,ijk 



(0) 
(1) _ 

p{uiUjUk + XijUk + XikUj + XjkUi - SijUk - SikUj - SjkUi), 



,(2) 



(9) 
(10) 

(11) 

(12) 



^Es.ijki = p[uiUjUiUk + (A^; - Sii)ujUk + (Ay - 5ij)uiUk + {\ik - Sik)ujUi 

+ {Xji - Sji)uiUk + {Xjk - Sjk)uiUi + {Xki - Ski)uiUj (13) 

+ i^ij ~ Sij){Xkl - Ski) + (Xik - Stk)iXji - Sji) + (A,;; - Sii){Xkj - 5kj)]- 



The body force term F{x, c, t) = g ■ V cf can be approximated as: 



TV 



^ (n-1) 

n=l ^ ' 



(14) 



Although Eq. ([2]) has an infinite Hermite expansion, a LB model is built under the assumption 
that to some level of approximation only the leading terms contribute explicitly to the hydro- 
dynamics (see |Shan et al. 20061 and thus the infinite series can be truncated. In this work, a 
4th-order Hermite expansion will be used. 

Given the above discussion, the ES-BGK equation Eq.Q can be rewritten in its truncated 
form, i.e., 



(15) 



For a 4th-order expansion, a discrete velocity set with 9th-order or higher accuracy is required 



Shan et al. 


2006 


). There are a number of ways of constructing a h 


Shan et al. 


2006 


Shan 


2010 


Chikatamarla & Karlin 


2009 


2006 



is to utilize the roots of the Hermite polynomial. In one dimension, the discrete velocities Cq. are 
just the roots of the Hermite polynomials, and the corresponding weights are determined by: 



,n~ll 



(16) 



Another useful procedure is to utilize the entropy construction (cf. Chikatamarla & Karlin 2006 



2009). Given one- dimensional velocity sets, those of the higher-dimension can be constructed 



using the "production" formulae (Shan et al. 2006). Once the discrete velocity set is chosen. 



the governing equation of the lattice ES-BGK model can be written as 

9/a 



^f<^ ^ (f f \ , 



(17) 



where /„ 



WgfEsjx.Cc.t) 



and go 



2.3 Remarks on the accuracy beyond Navier-Stokes 

As the distribution functions / and Jes are approximated and evaluated at the Gauss-Hermite 
quadrature points, cf Eq.([5]), we can relate the chosen Gauss-Hermite quadrature to the approx- 
imated moments of the distribution function. Generally speaking, more discrete velocities mean 
higher-order moments can be obtained more accurately. With sufficiently accurate quadrature 
and retaining increasingly high-order terms of the expanded fssy the LB model is becoming a 



better approximation of the original ES-BGK equation (cf Fig 20 1 
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When the Chapman-Enskog expansion is vaUd, we can evaluate model accuracy in terms of 
the Knudsen number. According to the Chapman-Enskog expansion, the distribution function 
/ is assumed to be given by an asymptotic series expanded in powers of the (formal) small 



parameter e (see Struchtrup 2005 & Chapter 4), 



(18) 



where e plays the role of the Knudsen number; f'^ is the Hermite approximation of the Maxwellian 



/m (its first four Hermite coefhcients can be found in Shan et al. (2006), cf. Eq.(3.12))). Fur- 

(19) 



thermore, the time and spatial variations are also measured in powers of e, i.e., 



and V = e V ( Shan et al. 2006 ) . Considering the form of / , f^g can be written as 



fN _ rCI 

Jes — J 



The shear stress aij should be also expanded as 



and can thereby be written as, 



where 



and 



fJV,(l) 
JES 



N,{2) 
ES 



UJ 



JES — J + <^JES + ^ JES ^ ' 
1, (1) (2) , 1/ (1) , (1) , (1) N (3) , 



,(2) 



(2), 



.(2), 



i3) 



(20) 



(21) 



(22) 



(23) 



(24) 



By matching the terms in the same powers of e, we can have the solution for the two leading 
orders, namely the NS order 



(0) 



and the Burnett order 



/ 



(2) 



JES 



(0) 



c-V+g.V,)/W+9i/°]- 



(25) 



(26) 



It is rather tedious to calculate the explicit form of /'^•' and Z'^'. However, here we only need 
to know their highest order as a Hermite polynomial of the particle velocity c. To determine the 
highest-order, we notice that the operators 5° and do not alter the order while the operators 
c-V and Vc increase the order by one. Therefore, the Burnett order solution f^^^^ is a polynomial 
of c of two-order higher than f^g. If desired, this discussion can be extended to solutions of 
higer order. Here we restrict the discussion to the first two orders. 

With the above discussion, we can estimate the accuracy level associated with the chosen 
expansion order and quadrature. As fourth-order terms of c are retained in the expansion of 
/eSj the Burnett-order solution of the distribution function can be represented by a sixth-order 
polynomial. For the heat flux, a third-order velocity moment, to be calculated with the Burnett 
level of accuracy, the chosen quadrature has to be accurate for an integration over the full space 



for a ninth-order polynomial. Hence, according to the Gaussian quadrature rule (cf., Shan et al. 



2006 Appendix in) , it is necessary to adopt a quadrature with an algebraic degree of precision 
beyond the ninth-order. However, as the ninth-order quadrature is not optimal for describing 



gas-wall interactions (Meng & Zhang 20116), an llth-order quadrature is chosen in this work. 



It is worth noting again that the above analysis relies on the validity of the Chapman-Enskog 
expansion. Consequently, one should be careful when interpreting the above conclusion for 
large Knudsen numbers. Furthermore, it may be better to understand the requirement on the 
quadrature for the corresponding accuracy as necessary rather than sufficient. Nevertheless, the 
above discussion is useful for understanding the capability of the LB model, which can be seen 
in numerical simulations below. 



6 



iJ+1 



* V * 



i+lj 



iJ-1 



y 



Figure 1: Schematic diagram of square lattices. 



3 Numerical implementation 

3.1 Scheme 

To solve Eq. 



(171 



various numerical schemes can be used, depending on the characteristics of 
the problem of interest. Typically, if a first-order upwind finite difference is chosen, we have a 
standard stream-collision scheme. 

As flows with various Knudsen numbers will be considered in the following simulations, the 
numerical scheme has to be carefully chosen to cope with discontinuities at the system boundaries 
caused by rarefaction (Hadjiconstantinou 2006). In particular, in order to intrinsically describe 



the gas-surface interaction, the kinetic boundary condition should be accurately implemented 
in the numerical scheme. Meanwhile, there is still no general and accurate implementation of 
the standard stream-collision scheme for high-order models, especially for rarefied flows. On the 
other hand, for high-order LB models, the discrete velocity points may not coincide with the 
lattice points, which makes the standard stream-collision procedure impossible. Therefore, in 
the interest of simplicity, a finite difference formulation has been chosen and more specifically, 
a Ist-order forward Euler method is coupled to a 2nd-order total variation diminishing (TVD) 
scheme ( |Kim et a/.[ |2008a| |SofoMa| |2009[ [Tbrol |2009[ ). 

The resulting scheme can be summarized as follows. Let /^f denote the distribution function 
value fa at the node (xi, yj) at the n-th time step (see FigjlJ. The distribution function update 
can be written as 



f 

•I a. 



f 

J a 



-^Q,i+l/2 -^0,1-1/2 



(27) 



j.n,j+l/2 _ T-n,i- 1/2 



y are the grid spacing in the x and y directions, respectively, and 



where 6^ and S. 

step; Cax and Cay denote the particle velocity components at the 
outgoing and incoming fluxes at the node (see Figjl]) are 



St is the time 
and y coordinates. The 
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(28) 

(29) 
(30) 
(31) 
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Figure 2: Schematic iUustration of waU boundary treatment. 



6" = 



•/ a, 2 


fn,3 
J a,i — l 
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J 


f"J _ 




j-nj + 1 

a, 2 


_ fn,j 
•taA 



(32) 
(33) 

(34) 



and the minmod flux Hmiter is given by 

* (6) = max[0,min(l,e)] 
For simphcity, only the formulae for Cay > and Cax > are given. 

3.2 Boundary condition 

Gas-wall interactions are captured by boundary conditions. The most popular boundary condi- 



tion is the Maxwell model (Cercignani 2000 Clerk Maxwell 1879), in which a fraction (1 — 7) 



of gas particles are assumed to undergo specular reflection while the remaining particles are 
diffusely reflected. In the limit of 7 = 1, the reflection becomes fully diffuse. 

As lattice models are derived from continuous kinetic equations, their boundary conditions 



may be obtained from the continuous kinetic theory ( Ansumali & Karlin 2002 1 . It has already 



been shown that lattice models with diffuse type wall boundary can describe a range of rarefied 



effects (see, e.g., Meng & Zhang 2011 & Kim e< aL 2008a Yudistiawan aZ. 2008 1. In this work, 
our primary interest is the model accuracy; as a result we will only implement the Maxwellian 
diffuse wall boundary. 



The boundary condition employed in this work is Version 1 of boundary conditions in Sofonea 



(2009), which is briefly described below. Firstly, the truncated Maxwellian distribution function 



can be written as 

/° = pS{T,u) 



+ 
+ 



pWa 

~6~ 
1 

24 

T- 1 



1 -I- CiUi + - [{ciUif - UiUt + {T ~ l)(ciQ - D)] 



(35) 



[(cjWi)^ - 3uiU, + 3{T - l)(ciCi - D - 2)] 

[(cjCi - D - 2)((wiCi)^ - UiU,) - 2{uiCi)'^] 
[(qc,)' - 2{D + 2)qc, + DiD + 2)]}. 
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As the discretization is conducted along a Cartesian coordinate system (see Fig[2|, the wall 
boundary condition can be written as 



fa,k — PW,kS{Tw,k, Uw,k) ' n > 0, 



(36) 
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Figure 3: The velocity and shear stress profiles for Couette flow of a hard-sphere gas at Ke, = 0.1, 
0.2 and 0.5. The symbols denote DSMC data and the lines represent the lattice ES-BGK results. 
The wall velocities ±0.1. 



where 

(Cc-n.)<0 

P^-f" ~ — Te ^1 Qfrr ~ \' v'^'-' 

2^ \ia- n\b{lw,k,uw.k) 

(€c-n)>0 

Here, subscript W denotes the computational nodes at the wall, pw,k is the density on the wall 
nodes k (see Figj2]), T^^^ is the temperature, the velocity and n the inward unit vector 

normal to the wall. The distribution function at the ghost nodes is assumed to be identical to 
those on the corresponding wall nodes. 



4 Validation 



We have validated the proposed scheme using a variety of flows, both in the low-Mach and high- 
Mach number limit. Unless otherwise stated, our results have been obtained using a llth-order 
of discrete velocity set, which represents a good compromise between accuracy (expected to be 
accurate at the Burnett level, see Sec 2.3) and computational efficiency. Our numerical results 



will be compared with DSMC and LVDSMC (see Homolle & Hadjiconstantinou 2007 Radtke 
& Hadjiconstantinou 2009 2011) simulations, while in some cases, numerical solutions of the 



ES-BGK equation and analytical results of the R13 model will also be presented. 

The flows considered here are one-dimcnsional and are confined by two infinite parallel plates 
in the x — z plane. The distance between the two plates (in the y direction) is L on which the 
Knudsen number 

Kd = ^J^Kn (38) 

is based. Comparisons with the hard-sphere gas will be performed by setting the Prandtl number 
to 2/3, i.e. b = —1/2. For the BGK gas, the Prandtl number is unity which corresponds to & = 0. 
Furthermore, the gas has the heat capacity Cp = 5/2, i.e., the spatial dimension D is set to be 
3, see Eqs. ^ and (35). Our simulation results {D = 3) are described below. 



4.1 Low Mach number flows 

Flows in the low-Mach number limit are of particular interest since they are most often en- 



countered in micro/nanofluidic applications (Cercignani 2006 Hadjiconstantinou 2006). Here 



we provide validation of our methodology using traditional Couette and Fourier flows. We also 
present solutions of an unsteady boundary heating problem that features coupled momentum 
and heat transfer and thus requires an accurate representation of the ratio of the associated 
transport coefficients, namely the Prandtl number. This problem also allows us to investigate 
the limitations of the proposed model in the presence of kinetic effects due to flow unsteadiness. 
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Figure 4: The temperature and heat flux profiles for a Fourier problem for a hard-sphere gas at 
Kd = 0.05, 0.1, 0.2 and 0.5 with the wall temperature difference 0.1. The symbols denote DSMC 
data and the lines represent the lattice ES-BGK results. 



4.1.1 Steady Couette and Fourier Flows 

Figures |3] and |4] show comparison between our lattice results and DSMC simulations for a Couette 
and a Fourier heat tranfer problem respectively. In the Couette flows, the walls move with 
velocities = ±0.1, while in the Fourier problems the wall temperatures are r„, = 1 ± 0.05. 
Fig. [4] shows that the proposed lattice ES-BGK model can be used to approximate the hard- 
sphere gas thermal conductivity by setting the Prandtl number to 2/3. Both figures show that 
the agreement in the slip-flow regime is excellent; as Kd increases into the transition regime, a 
small amout of error is evident. This error increases further as Kj^ increases, but remains very 
reasonable for Kd < 0.5; for example, for both flows, at Kd = 0.5 the error is less than 5%. 



4.2 Unsteady boundary heating problems 

In this section we present results from an unsteady boundary heating problem of the kind studied 



in (Manela &: Hadjiconstantinou 2007 1 ; these studies are motivated by the common occurrence 
of applications involving time-varying boundary temperatures in micro-electro-mechanical de- 
vices. In this problem, the two walls confining the gas are heated uniformly via prescribed time 
dependent temperatures of the form = 1 + F{t). Here, F{t) is taken to be a sinusoidal func- 
tion 8m{9t) with an amplitude (0.002 in simulations) that is sufficiently small to justify a linear 
assumption. As this is a time-varying problem, another non-dimensional parameter needs to be 



introduced (Manela & Hadjiconstantinou 2008 2010 ), namely the Strouhal number, defined as 



St = 9- 



0L 



(39) 



Kinetic effects become important as both the Strouhal and Knudscn numbers increase as shown 



in Fig. 6 in Manela & Hadjiconstantinou (2010) 



Here, simulations will be conducted for a range of Strouhal and Knudsen numbers. Prandtl 
numbers corresponding both to the hard-sphere and BGK model will be considered. The accuracy 
of the lattice results will be evaluated by comparison with simulation results obtained using the 



recently developed low- variance deviational simulation Monte Carlo (LVDSMC) method (Ho- 



molle & Hadjiconstantinou 2007 Radtke & Hadjiconstantinou 2009[ 20111 because simulation 



of this low-Mach number ffow by DSMC is prohibitively expensive. The LVDSMC method was 
developed to specifically address this DSMC limitation, namely the prohibitive computational 
cost associated with low-Ma and more generally low-signal flows. LVDSMC provides efficient 
numerical solutions of the Boltzmann equation for all Knudsen numbers ([Wagner 



20081 and 



arbitrarily small signals by simulating only the deviation from equilibrium. This control- variate 
variance-reduction formulation introduces deterministic knowledge of a nearby equilibrium state, 
and thus significantly reduces the stastistical uncertainty associated with the Monte Carlo ap- 
proach without introducing any approximation. The resulting computational benefits in the limit 
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Figure 5: The velocity and temperature perturbations for a hard-sphere gas subject to a sinu- 
soidal heating at i = 37r/2 at different Kn and St. The symbols correspond to the data of the 
lattice model (Pr — 2/3) and the Hues are the results of the LVDSMC method. 




Figure 6: The velocity and temperature perturbations (AT) for a hardsphere gas subject to a 
sinusoidal heating a.t t — Sir/ 2 at different Kn and St. The symbols correspond to the data of 
the lattice model {Pr = 2/3) and the lines are the results of the LVDSMC method. 



of low speed flows such as the one considered here are typically very large (Homolle & Hadji- 



constantinou 



20071 



Following Manela & Hadjiconstantinou ( |2010[ ), comparison with LVDSMC results will be 

(40) 



performed using the following two Knudsen numbers: 

16 



5V27r 



Kn, 



for the hard-sphere gas with Pr = 2/3, and 



Kb = 



Kn 



(41) 



for the BGK gas with Pr = 1. Details on the LVDSMC simulations of the problem studied here 



can be found in Manela & Hadjiconstantinou (2008 2010). 



The results for various Strouhal and Knudsen numbers for the hard-sphere and the presented 
in Figsj5|8j where both the velocity Uy and temperature perturbation AT — T—1 are normalized 
by the amplitude of F{t). Figs 9p2 show results for the BGK model. Symbols denote the lattice 
model proposed here, whereas lines denote LVDSMC data. At St = tt\/2/16, the two sets of 
results show good agreement, even when Ke = 0.5 or Kb — 0.5. With increasing Strouhal 
number, disprepancies between the two appear and become larger. For St = 7r-\/2/4, significant 
disagreement is observed even for Ke as low as 0.05. Larger Strouhal numbers lead to stronger 
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Figure 7: The velocity and temperature perturbations for a hard-sphere gas subject to a sinu- 
soidal heating at i = Sir/ 2 at different Kn and St. The symbols correspond to the data of the 
lattice model {Pr = 2/3) and the lines are the results of the LVDSMC method. 
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Figure 8: The velocity and temperature perturbations for a hard-sphere gas subject to a sinu- 
soidal heating at i = Sir/ 2 at different Kn and St. The symbols correspond to the data of the 
lattice model {Pr = 2/3) and the lines are the results of the LVDSMC method. 




Figure 9: The velocity and temperature perturbations for a BGK gas subject to a sinusoidal 
heating at t = 37r/2 at different Kn and St. The symbols correspond to the data of the lattice 
model {Pr ~ 1) and the lines are the results of the LVDSMC method. 
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Figure 10: The velocity and temperature perturbations (AT) for a BGK gas subject to a sinu- 
soidal heating a,t t — 3tt/2 at different Kn and St. The symbols correspond to the data of the 
lattice model {Pr ~ 1) and the lines are the results of the LVDSMC method. 




Figure 11: The velocity and temperature perturbations for a BGK gas subject to a sinusoidal 
heating at t = 37r/2 at different Kn and St. The symbols correspond to the data of the lattice 
model {Pr = 1) and the lines are the results of the LVDSMC method. 




Figure 12: The velocity and temperature perturbations for a BGK gas subject to a sinusoidal 
heating at t = 37r/2 at different Kn and St. The symbols correspond to the data of the lattice 
model {Pr ~ 1) and the lines are the results of the LVDSMC methods. 
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Figure 14: The velocity and temperature profiles (AT = T — T^) for Couette flow of a hard- 
sphere gas at Ko = 0.1, 0.2 and 0.5. The symbols denote the DSMC data and the red lines 
represent the lattice ES-BGK results. The wall velocities are Uy^ = ±0.2. For further comparison, 



the temperature profiles predicted by the R13 model (Taheri et al. 2009) are also presented with 



the black lines where the line styles same to those of the lattice ES-BGK results are used to 
distinguish the Knudsen numbers. 



rarefaction effects, so this disagreement can be attributed to the moderate discrete velocity 
set. Overall, with this moderate discrete velocity set, the present model can give reasonable 
predictions for flows with a Knudsen number up to 0.5 and a Strouhal number up to 7r-\/2/8. 
If highly accurate results are desirable, more discrete velocities are required, leading to higher 
computational costs. 

4.3 High-Mach number flows 

In this section we present simulation results for high Mach number flows such as Couette flows 
with Uw ± 0.2 and forced Poiseuille flows with non-dimensional forcing magnitude g = 0.22. 
Although in the Couette flows the wall speed is only a factor of 2 larger than the flows examined 
in section |4.1[ our focus here turns to the resulting temperature field that, despite the small 
temperature differences involved, reveals useful information. 

Fig. [13] shows a comparison between our lattice model, DSMC and numerical solution of 
the ES-BGK equation for a Couette fiow at Kd = 0.05. The comparison reveals that at these 
small Ku, the lattice model can capture both kinetic (e.g. slip/temperature jump) and non- 
equilibrium effects quite accurately. For comparison, the results of the BGK equation and its 
lattice version (6 = 0) are also presented. These two models predict a temperature maximum 
that is 30% higher than the hard-sphere result. 
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Figure 15: The shear stress {(Txy) and heat flux (qy) profiles for Couette flow of a hard-sphere 
gas at Kjj = 0.1, 0.2 and 0.5. The symbols denote DSMC data and the lines represent the lattice 
ES-BGK results. The wall velocities are = ±0.2. 
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Figure 16: The shear stress {cyy) and heat flux [q^) profiles for Couette flow of a hard-sphere 
gas at Ko = 0.1, 0.2 and 0.5. The symbols denote direct solution of the ES-BGK equation, the 
red lines represent the lattice ES-BGK results and the black lines are those of the R13 model 
(Taheri et al. 2009). The Knudsen numbers for the R13 model are denoted using the same line 
styles as the lattice ES-BGK results. The wall velocities are u^, = ±0.2. 




Figure 17: The velocity and temperature profiles (AT = T — T^) for the force-driven Poiseuille 
flow of a hard-sphere gas at Kd = 0.1, 0.2, 0.4 and 0.5. The symbols denote DSMC data and 
the lines represent the lattice ES-BGK results. 
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Figure 18: The shear stress and heat flux profiles for the force-driven Poiseuille flow of a hard- 
sphere gas at Kd — 0.1, 0.2, 0.4 and 0.5. The symbols denote DSMC data and the lines represent 
the lattice ES-BGK results. 
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Figure 19: Comparisons of the velocity and temperature profiles (AT — T — T^) for the force- 
driven flow among the ES-BGK equation, lattice ES-BGK model and DSMC method at Ku = 0.1 
and 0.5. 




Figure 20: Convergence comparison for a Couette flow at Kjj = 0.2 with wall velocities = 
±0.2. Lines denote the lattice model data, and the adopted discrete velocity order as denoted in 
the legend. Improved accuracy is observed as the number of discrete velocities is increased. 
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Figure 21: Comparisons of the velocity and temperature profiles for the combined Couette- 
Fourier flows. The wall velocities are = ±0.2 while the wall temperature difference is 0.3. 
The lines are the lattice ES-BGK model data and the symbols are the results of the ES-BGK 
equation. 



Results for Coucttc flows in the transition regime are presented in Figs. 14 and 15 Overall, 



our model predictions for the velocity, shear stress and heat flux are close to those of DSMC 
even for Kj^) = 0.5. However, our predictions of the viscous-heating-induced temperature fleld 
start to deviate from the DSMC results at Kd = 0.2. Remarkably, the heat flux predictions of 
the lattice model are still in excellent agreement with the DSMC data for Ku as large as 0.5. 

To further validate the model accuracy beyond the NS level, the profiles of shear stress ayy 
and stream- wise heat flux are presented in Fig. (16), which are both zero at the NS level 
of approximation. However, our lattice model can give reasonable predictions for ayy up to 
Kd = 0.2 which is consistent with the temperature prediction. For the heat flux q^, we have 
similar observation, which further conflrms that our model can describe kinetic effects beyond 
the NS level with the chosen moderate 11-th order quadrature, which is consistent with the 
discussion in Sec|2.3| It is also interesting to compare with the predictions of the R13 model. 



which is supposed to give a stable set of transport equations of the super-Burnett order ( Taheri 



et al. 2009). Fig. (16 1 shows that the R13 model and our lattice model give similar results. 



The simulation results for Poiseuille flows at four Knudsen numbers are depicted in Figs. [17] 
and 18 This flow is more numerically challenging compared to the Couette flow. Even for Kd as 
small as 0.1, the temperature proflles show signiflcant difference between the lattice model and 
DSMC results, even though the bimodal temperature distribution (Mansour et al. 1997) can 
be qualitatively captured by the lattice model at Kjj — 0.2. As the Knudsen number increases 
further, the difference becomes even larger, though the temperature minimum at the center of 
the computational domain can be predicted qualitatively. Remarkably, but also consistently with 
the Couette flow case, the heat flux proflles show good agreement with the DSMC data despite 
the more significant discrepancies in the temperature profiles. 

The inability of our model to predict accurate viscous heating may be caused by two fac- 
tors. First, as shown in the derivation process, the lattice model approaches the corresponding 
kinetic model via the moment representation with increasing order of Gauss-Hermite quadra- 
ture. However, as discussed in Sec 2.3 the chosen moderate discrete velocity set is only expected 
to approach the Burnett level, which may not be sufficient to capture all kinetic effects in the 
presence of a non-negligible Mach number. This is verified in Fig. [20] where its is shown that 
enriching the discrete velocity set can improve the model accuracy. 

In addition to the error due to a finite number of discrete velocities, the ES-BGK equation 
itself is a model which may be failing in these flows. In fact, as shown by some studies (e.g 



Gallis 



& Torczynski |2011J , the ES-BGK equation tends to produce inaccurate velocity distribution 



functions. The possibility that the ES-BGK model is itself contributing to the error can be 
confirmed by the results presented in Fig. [19] At = 0.5 the lattice model result shows 
a significant deviation from the numerical solution of the ES-BGK model. Meanwhile, the 
temperature profile from numerical solution of the ES-BGK equation itself shows a large deviation 
from the DSMC result. 
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While the comparisons with the DSMC solutions for the temperature difference of 0.1 have 
been shown in Fig. |4| a larger temperature difference of 0.3 is investigated for the combined 
Couette-Fourier flow. The results are shown in Fig. 21 where excellent agreement is observed 
between the lattice ES-BGK model and the ES-BGK equation. 



5 Concluding remarks 

We have presented and validated a lattice Boltzmann model using a systematic Hermite mo- 
ment representation of the ES-BGK equation. The resulting lattice ES-BGK model features an 
adjustable Prandtl number and may thus be more appropriate for use in coupled thermofluidic 
phenomena. This generic procedure may be applied to other kinetic model equations. 

We have validated the lattice model for combined thermal and rarefaction effects by comparing 
its predictions with the DSMC and LVDSMC results as well as numerical solutions of the ES-BGK 
model. We find that, for a moderate llth-order discrete velocity set, the proposed model can 
provide reasonable predictions for Couette and force-driven Poiseuille flows for Knudsen numbers 
up to 0.5. In addition, it is able to accurately predict heat conduction in the slip-flow and early 
transition regimes. This finding extends to unsteady problems provided the additional kinetic 
effects due to their time-dependence are also considered: for an unsteady boundary heating 
problem reasonable agreement with LVDSMC simulations is observed for Knudsen numbers up 
to 0.5 and Strouhal numbers up to 7r-\/2/8. 

The solutions obtained by the lattice Boltzmann method are expected to approach the "true" 
solutions of the ES-BGK equation as the number of discrete velocities and the Hermite expansion 
order is increased. However, the moderate discrete velocity set used here already represents a 
reasonable compromise between computational efficiency and modeling accuracy for flows with 
a range of Knudsen and Strouhal numbers. Although we have tested the model for rarefied gas 
thermal flows, this model, in principle, can be used for liquid thermal flows, which will be the 
subject of future work. 
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